print("########################################################################################")
print("# Title: Sub-network GSA                                                               #")
print("# Author: Yongkee Cho (yongkeecho@wustl.edu)                                           #")
print("# Date: Aug. 27 2012                                                                   #")
print("# Description: This script runs probeset-level SAM analysis and                        #")
print("#              returns gene-level SAM-t score (best among matching                     #")
print("#              probesets                                                               #")
print("# Requirements: GSA, bioC/org.Hs.eg.db                                                 #")
print("#                                                                                      #")
print("# Usage: Rscript subnet_gsa.r test.gmt test_xg.csv 'Two class unpaired' 0.5 output.gmt #")
print("########################################################################################")

# reads arguments
args <- commandArgs(trailingOnly = TRUE)

# mapping candidate gene symbol to ENTREZ gene ids
print("Reading sam data file...")
if ("org.Hs.eg.db" %in% rownames(installed.packages()) == FALSE) {
   source("http://bioconductor.org/biocLite.R")
   biocLite("org.Hs.eg.db")
}
library(org.Hs.eg.db)
x <- org.Hs.egSYMBOL
mapped_genes <- mappedkeys(x)
xx <- as.list(x[mapped_genes])

# GSA analysis
if ("GSA" %in% rownames(installed.packages()) == FALSE) {
   install.packages("GSA")
}
library(GSA)
subnets <- GSA.read.gmt(args[1])
xg <- read.table(args[2],header=TRUE,sep=",",quote="\"")
geneidg <- rownames(xg)
y <- as.integer(substr(colnames(xg),2,2))
gsa <- GSA(xg, y, genenames = geneidg, genesets = subnets$genesets, minsize = 3, maxsize = 100, resp.type = args[3], nperms = 100)
nn <- length(subnets$genesets)
subnets_names <- rep("", nn)
for (i in 1 : nn) {
  subnets_names[i] <- toString(xx[subnets$genesets[[i]]])
}
fdr <- as.numeric(args[4])
gsa_list <- GSA.listsets(gsa, geneset.names = subnets_names, FDRcut=fdr)
sig_subnetidx <- as.integer(c(gsa_list$positive[,1], gsa_list$negative[,1]))
sig_gmt <- readLines(args[1])[sig_subnetidx]
writeLines(sig_gmt, args[5])
